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Probing magnetic fields in volume 

with multi-frequency polarized synchrotron emission. 
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ABSTRACT 

We investigate the problem of probing the local spatial structure of the magnetic field of 
the interstellar medium using multi-frequency polarized maps of the synchrotron emission 
at radio wavelengths. We focus in this paper on the three-dimensional reconstruction of the 
largest scales of the magnetic field, relying on the internal depolarization (due to differential 
Faraday rotation) of the emitting medium as a function of electromagnetic frequency. We 
argue that multi-band spectroscopy in the radio wavelengths, developed in the context of high- 
redshift extragalactic HI lines, can be a very useful probe of the 3D magnetic field structure 
of our Galaxy when combined with a Maximum A Posteriori reconstruction technique. 

When starting from a fair approximation of the magnetic field, we are able to recover 
the true one by using a linearized version of the corresponding inverse problem. The spectral 
analysis of this problem allows us to specify the best sampling strategy in electromagnetic fre- 
quency and predicts a spatially anisotropic distribution of posterior errors. The reconstruction 
method is illustrated for reference fields extracted from realistic magneto-hydrodynamical 
simulations. 



1 INTRODUCTION 



The problem of studying the magnetic field structure of our 
Galaxy using measurements of the synchrotron emission of high 
energy electrons in the Galactic magnetic field is an old one 
dOinzburg & SvrovatskiillT965l : iRuzmaikin etalJI 19881 : [Beck et al.l 
Il996h . The fact that the emitting medium is itself magnetized 
induces a differential Faraday rotation of the different emission 
planes transverse to the line of sight, resulting in a well known de- 
polarization effect of the integrated emission that depends strongly 
on the el ectromagnet ic frequency. This effect, described in the first 
. !^ ' place bv lBurnl ( ll966h in the case of a constant magnetic field, has 
, been further studied in semi-analytically for given functional forms 
5^ of the magnetic field; it has also been studied fro m the statistical 
, point of view in some asymptotic regimes (see e.g. ISokoloff et al.l 
Il998l) . In the present work, we want to consider the more am- 
bitious problem of using this depolarization effect, together with 
the solenoidal character of the magnetic field, to recon.struct the 
magnetic field structure from a set of polarized maps of the syn- 
chrotron emission of an ionized medium at different electromag- 
netic frequencies. With the upcoming prospect of detailed Multi- 
band spectroscopy in the radio wavelengths jRottgerind bool 
[purlanetto &. Briass 2004), developed in the context of Galactic 
and high-redshift extragalactic HI lines, this type of investigation 
should become possible. 

A statistical inference of the measurement of the Galac- 
tic magnetic field correlator as a function of scale from multi- 
frequency polariz ation measurements h as already been success- 
fully achieved by IVogt & Enfilii] { 2005h in the case of the Fara- 
day rotation of the polarized light from background objects by the 



intra-cluster magnetized plasma. In this case, there is no depolar- 
ization effect due to differential Faraday rotation, and the relation- 
ship between the measured polarization at a given frequency and 
the polarization of light in the source plane is linear in the (longitu- 
dinal) magnetic field strength. The linearity of the problem makes 
the statistical analysis tractable in the former case. In the case that 
we investigate, the emitting and the rotating medium are the same, 
which results in depolarization effects of the emitted light. More- 
over, the synchrotron emissivity itself depends non-linearly on the 
field strength transverse to the line of sight. The reconstruction of 
the magnetic field structure from the polarization data is in this case 
a non-linear inverse problem. Finally, we must note that to address 
the full problem of reconstruction of the magnetic field from the 
depolarized synchrotron emission we need in principle knowledge 
of both the thermal electron spatial distribution Uc and the spa- 
tial distribution of cosmic ray electrons nr, when, in comparison, 
the inference of the magnetic energy spectrum from the rotation 
measures of background sources only requires knowledge of the 
thermal electron distribution. 

In a first attempt at reconstructing the magnetic field, and for 
the sake of clarity, we make the assumption that the fluctuations of 
the thermal and cosmic ray electrons can be neglected compared 
to the fluctuations in the magnetic field itself. This assumption, if 
physically unrealistic, allows us to show the specific influence of 
the magnetic field statistical properties on the quality of the re- 
construction. In the first sections, we thus consider the electronic 
distributions (both thermal and relativistic) as constant, and discuss 
the reconstruction of the magnetic field using only the leading cou- 
pling coefficient in the equation of radiative transfer. In the (thin 
medium, strong rotativity) limit that we assume for this work, this 
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leading term is tlie usual Faraday term, responsible for the rotation 
of the plane of polarization. We will assume that the Faraday co- 
efficient is dominated by the thermal electrons, which is a reason- 
able assumption in non-relativistic astrophysical plasmas. Finally, 
in section|4l we relax the unrealistic assumption of a constant ther- 
mal electrons density, and show that our method can still be used to 
reconstructed the magnetic field when the electronic density is spa- 
tially varying but known a priori, using simulated data sets from a 
magneto-hydrodynamical (MHD) simulation. 

This paper is organized as follows: in section |2] we discuss 
the fonctional dependence of the polarization of the synchrotron 
emission and its variation with electro-magnetic frequency on the 
underlying magnetic field. We present a discretized version of this 
functional dependence that will be useful in the context of the re- 
construction from discrete polarization data. In section[3]we inves- 
tigate the reconstruction of the magnetic field from simulated multi- 
frequency polarized data, when the functional dependence on the 
magnetic field has been linearized around a "mean" field. Taking 
advantage of the linear nature of this approximate problem, we give 
a strategy for choosing the best electromagnetic frequencies of ob- 
servation, and investigate the statistical anisotropy of the magnetic 
field reconstruction errors. Finally, in section U we investigate the 
validity of the linearization procedure used in the precedent section, 
as a function of the quality of our prior knowledge of the magnetic 
field structure. We show how the approximate, linearized inverse 
problem investigated in this work could be used as a building block 
of the fully non-linear reconstruction problem. We emphasize that 
any gradient-based non-linear minimization algorithm can be de- 
composed into linear sub-problems, thus justifying the study of the 
linearized problem. In this context, we investigate how the condi- 
tioning of the linearized problem varies with the properties of the 
reference magnetic field around which the problem is being lin- 
earized. In particular it is illustrated on a realistic reference field 
from a MHD simulation. Finally, using the same MHD simulation 
data, we show that our method can deal with a non-constant elec- 
tronic density, provided it is known a priori. In section|5] we sum- 
marize the main results of the paper, recalling the main simplifying 
assumptions used to derive them (notably the assumed-known elec- 
tronic density hypothesis) and discuss how this assumption could 
be possibly alleviated by additional data (e.g., H^, free-free) or by 
using second-order coupling terms involving the circular polariza- 
tion in the case of relativistic sources (see[C}. We conclude on how 
the different results of the paper could be used to tackle the fully 
non-linear reconstruction of the magnetic field. 



2 POLARIZED EMISSION 

Our objective is to recover the magnetic field given observed po- 
larization maps at different wavelengths. We tac kle this ill-posed 
problem by means of an inverse problem approach jTarantolal 19871) 
which involves recovering the magnetic field that gives a polariza- 
tion consistent with the observations while obeying some a priori 
properties. These priors are strict constraints, such as V ■ B = 0, 
to insure that the sought field is physically meaningful and a reg- 
ularization to lever the degeneracies of the inverse problem while 
avoiding artifacts due to noise amplification. We first derive the di- 
rect model of the polarization given the magnetic field and then 
introduce the inverse problem approach in a Bayesian framework. 



2.1 Direct model 

We only consider here the Faraday rotation in the transfer equation, 
and neglect all other coupling terms. In this case, the transfer equa- 
tion of the Stokes parameters of linear polarization (Q, U) can be 
integrated formally. We assume here that the density of electrons 
is constant, or that its fluctuations are only important on scales that 
are not considered here. 

Consider a slab of ionized magnetized medium of width L 
which is emitting synchrotron radiation. The polarized emission, 
as a f unction of frequency , integrated over the line of sight then 
reads jSokoloff etaljl998l) : 



P = Q + W ^ / e(r) 



2iV(r) 



Az , 



(1) 



with Q and U are the usual Stokes parameters, e(r) the synchronton 
emissivity which obeys: 



e(r) = An,{r) |Bx(r)|^ //""^ 



(2) 



and ^{r) the sum of the Faraday rotation and the primordial orien- 
tation: 



ip(r) ^ tt/2 + aictan (By/ Bx) A ^ / UcB^dz, 



(3) 



where r = (a;, y, z) = (xx, z) is the coordinate in the slab, v is 
the frequency, and B = {Bx, By, B^) — (Bx, Bz) is the magnetic 
field. In equation K reads: 

3 

K = , (4) 



while, in equation l|2j, A is given by 

7-1 



A = 



3(?e 



16 TV eo rrio c \2-Km^c* J ^\ 12 12 

where Eo is the energy scale of the relativistic electron spectrum, 
iric and go stand for the mass and the charge of the electron, ric 
and fir are the thermal and relativistic electron densities supposed 
constant, while the exponent 7 stands for the spectral index of the 
cosmic ray electrons, c is the speed of light, eo is the electric per- 
mittivity and F is the Euler gamma function. The lengths are in 
kilo-parsec (kpc) and so the density in kpc~^, the magnetic fields 
in micro-Gauss (fiG) and the frequencies in giga-Hertz (GHz). Re- 
expressing the intrinsic polarization phase in terms of powers of the 
magnetic field components, we get the following expression for the 
polarization: 

_i -7—3 
P(xx,!/) = Aj/-"^ / n,(xx,z) (B^ + B^)~ (xx,2:) 

J —00 

X (B^-B^+2iB,B,,)(xx,^) 

X exp^^^^ J (noB2)(xx, z") dz"^ dz . (5) 

As real data come in discrete form, let us discretize this expression 
by replacing all integrals with sums, assuming a regular discretiza- 
tion grid that will be defined more precisely below. Equation Q 
then reads 

P(xx,!/) ^^/ii^""^^ n,(xx,z) (B^ + B^)^ (xx,2) 



X (B^-B^-f 2iB,B„)(xx,z) 
'2iKh 



X exp 



Quiz' - z) {Ua B,)(XX, z'] 



(6) 
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Here Oh is the Heaviside function (6h{x) = 1 for a; > and 
elsewhere), and h the discretization length along z. Equation ^ is 
formally a function of B = {{Bx, By, -Bz)}^ where we use bold 
symbols to represent the discretized vector fields and r is a triple 
index spanning the magnetized volume on a regular cubic mesh 
with cell size h. 

The solution to the inverse problem will be obtained by means 
of minimization of some merit function (as explained in what fol- 
lows), we therefore need to compute the partial derivatives of the 
polarization with respect to the magnetic field. Let us first com- 
pute the derivatives with respect to the transverse components of 
the field: 



dBx{v') 
'1 



^5,,{v-v')An,{v')hv-'^{Bl+Bl)^{v') 



^-Bl + l- 



BlBx + -i{-i-l)BlBy + 2\B. 



X exp [ ^^^^ O-a {z" - z) {n^ B^ ) (x^ , z'' 



(7) 



with r = (xx,z), r' — (x'^jz') and So Dirac's delta function. 
The derivative with respect to By follows closely, with the square 
bracket term becoming: 



1+7 

^By 



7- 



^Bl B,+i(7-l)B, Bl + 2iBl 



(8) 



which corresponds to a 7r/2 rotation in the plane perpendicular to 
the LOS. We see that in both cases the phase term is unaffected 
since it is only a function of the longitudinal magnetic field com- 
ponent B^. Finally let us compute the derivative with respect to 
S,: 

'-^^ = &(xx-xl)2iX^/.^."^ 



'2iKh 



{Bl-^ B^y + 2iBxBy)ix'^,z) 



x9h{z' - z) exp — — 5Z ^h(2" - 2)(ne-Bz)(xi 



(9) 



We note that here the phase term, not the emissivity layer term, is 
involved. The case 7 = 3 is detailed in Appendix IaI and leads to a 
simplification of the above equations. 

2.2 Maximum A Posteriori formulation 

From the direct model, we can express the observed data as: 

d™^ = P((xx,!/)m,B) +e™^, (10) 

with m an index which spans the mixed frequency position-on-the- 
sky cube, (xx,i^)m the corresponding coordinates, B the actual 
magnetic field and em an error term which accounts for noise and 
model approximations. Using vector notation, equation dlOb simpli- 
fies to: d = P(B) +e with d = {dm} the vector collecting all the 
observations, P(B) = {P((xx, !^)m, B)} and e = {em}- Our 
inverse problem is to recover the magnetic field vector, B, given 
some noisy measurements of the polarization, d. Due to the un- 
known errors in equation dlOt and to possible strict degeneracies 
of the direct model, there is not a unique magnetic field that yields 
a polarization consistent with the observations. We therefore need 
some means to select a unique solution and, hopefully, the best one 
given the data. 

Probabilities provide a consistent framework to define such 



a solution; we thus define the sought magnetic field as being the 
most likely given the observations. It is the one which maximizes 
the posterior probability: 



B^ 



argmaxP(B|d) 

B 



(11) 



and whi ch is termed as the maxim um a posteriori (MAP) solution 
(see e.g. lPichon & Thiebaull 19981) . By Bayes' theorem, P(B|d) = 
■p(d|B)7'(B)/7'(d), and since ^(d) does not depend on the 
sought parameters B, this amounts to maximizing ^(djB) 'P(B). 
The term P(d|B) is the likelihood of the data given the model, 
while the term ^(B) accounts for any a priori knowledge about 
the magnetic field. We can anticipate two types of priors: (i) the 
strict constraint that, to be physically meaningful, the field should 
be solenoidal: VB = 0; (ii) some so-called regularization con- 
straint to overcome the ill-conditioning of the inverse problem and 
to enforce the unicity of the solution. Without loss of generality, we 
state that the probabilities writes: 

V{d\B) = Ki exp(-i £(B)) , (12) 
K2 exp(-i7^(B;/x)) , if VB = 0, 



P(B) 







otherwise. 



(13) 



where the factors ki and K2 do not depend on B and /i accounts 
for parameters to tune the regularization. Finally, taking the log- 
probabilities and discarding constants, the maximum a posteriori 
magnetic field writes: 



with: 



Bmap = argmin Q(B) , 

B,VB=0 



Q(B) = £(B) + 7^(B;/i). 



(14) 



(15) 



which is the objective function. Before going into the details of the 
expressions of £(B) and 72.(B;/i) we can already note that the 
solution Bmap will depend on the data d and on the regulariza- 
tion parameters ^. The value of ^ can be chosen, e.g., to provide 
the best bias-variance co mpromise on the sought solution jWahbal 
ll99(]l : lGolubetalJ200d) . 



2.2.1 Likelihood 

Assuming Gaussian statistics for the noise and model errors, the 
likelihood of the data is the so-called and writes: 

£(B) = (d-P(B))^-C-^-(d-P(B)) (16) 

with Cn the covariance matrix of the errors. There is a slight issue 
here because we are dealing with complex values. Since complex 
numbers are just pairs of reals, complex valued vectors such as d, 
P(B) and e can flattened into ordinary real vectors (with dou- 
bled size) to use standard linear algebra notation. This is what is 
assumed in equation l ll6t . Under these conventions, the covariance 
matrix of the errors writes Cn = (e ■ e^) with ^ to denote trans- 
position. 



2.2.2 Regularization 

The regularization term 7?.(B; /i) implements loose constraints to 
avoid over-fitting the data and enforce local unicity of the solution 
(see section |43] |. Requiring that the magnetic field be as smooth 
as possible (while being consistent with the data) matches these 
requirements and is supported by physics since the magnetic field 
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should have no discontinuities. To simplify further computations, 
we choose the following particular expression of the regularization 
TZ to favor the smoothness of the field: 



7^ = ^i.||A' 



«/4, 



OC 



E 



kriBi- 



(17) 



which scales as the integrated norm of the spatial Laplacian of the 
field to the power q/4. For a periodic field, this generic smoothing 
penalty is diagonal in Fourier space. In addition, if the model B is 
Gaussian and scale invariant, then a may be chosen to be the power 
law index of the power spectrum |B|| of the field. In this case, 
choosing the specific value of the hyperparameter, /is = l/|B|^^i, 
the MAP solution correspond to the minimal variance Wiener fil- 
tered data. 



2.2.3 Imposing VB = 

For simplicity, we assume here that the magnetic field is multi- 
periodic, with period L in all three directions. We may then rewrite 
the magnetic field as: 

-1 



B = F 



(Bxiexi + 5x26x2) = n-B, 



(18) 



where = /N^ and F is the forward DFT operator, (e|| = 
k/|k|, exi, 6x2) form a spherical basis in Fourier space, while 
-B±,i, i=l,2 are the projections over that basis of the Fourier com- 
ponent, B = F ■ B of the field. Equation J18t defines the projector 
n — F~^ • (ex ® ex) • F. Such a field satisfies by construction 



k ■ B = , which implies V • B = . 



(19) 



In fact, there is a slight complication at the Nyquist frequencies 
where only one component of the field is free, see appendix iBl 

Note that the diver gence free condition could also be imposed 
by other means (see e.g. lNocedal & Wrighll2006l) . For instance, by 
adding a quadratic penalty term like ^^{^ B)^ to the total penalty 
Q(B). We however found that, in practice, the projector II led to 
a better conditioned reconstruction problem. 



2.3 Implementation 

Given equations l ll6b and l ll7b the objective function writes: 
Q = (P - d)t . C-^ ■ (P - d) + II A"/^Bf . 



(20) 



To minimize Q(B), we used a variable metric limited m emory 
optimization method with BFGS updates (Nocedal' 1980') called 
VMLM and implemented in OptimPaclQ ( Thiebaut 2002). Finding 
the optimal solution, equation ( 114b . involves computing the gradi- 
ent of equation ( 120b with respect to B. Now differentiating equa- 
tion ( 116b with respect to a magnetic field components we get 



OB, 



= 27^e 



(21) 



where dP/dBi for i — x,y,z are given by equations ([7]l and ([9]l. 
Similarly, differentiating equation ( 117b with respect to B yields 

dn 



^ = MsF-^.B.|kr 



(22) 



The VMLM algorithm is a quasi-Newton method which proceeds 
by solving successive linear problems. Let us therefore first con- 
sider in the next section a linearized version of our inverse problem, 

^ OptimPack is freely available at http://www-obs.univ- 
lyonl . fr/labo/perso/eric . thiebaut /opt impack . html. 



which may correspond to a physically motivated problem when a 
good first guess for the magnetic field is known. 

Note finally that equations (Q and ^ imply that dx^ /dBi — 
at Bx = By = 0. Note also that if {Bx, By, B^) is a. solution 
to equation ([5j, so is {—Bx,—By,Bz). Consequently we expect 
that the will be strongly multivalued as a function of B0. The 
smoothing penalty should in part prevent a pixel-by-pixel flip of 
the X and y component. It remains nonetheless to be shown that the 
zero divergence condition is sufficient to avoid flipping the field in 
regions bound by zeros of these two components, if such regions 
exist. Addressing these issues will be the topic of another paper. 



3 LINEARIZATION 

Let us first consider the situation when a fairly good guess for 
the overall magnetic field. Bo, is known, on the basis, say of 
a first large scale investigation, or via some modelling of the 
field as a function of t he underlying density (e.g. ICao et al]|2006l : 
iKachelrieB et alj|2007h . Let us then seek the departure from this 
guess. It is then legitimate to assume B — Bo + SB, with, possi- 
bly (if the prime guess is accurate enough) <5B/jBo| <^ 1, so that 
equation (|5) becomes: 



5P : 



(23) 



where the tensor dP/dBi is given by its components, equations 
(0, (ID and while SP = P - P(Bo). Now equation ^ is 
likely to be a much better behaved equation as the linearity warrants 
convexity of the objective function, hence the formal unicity of the 
solution. 

In this paper, we will address two linear problems in turn, one 
of academic interest, to understand the properties of the inverse 
problem at hand, while the second one should allow us to carry 
realistic reconstructions, in the regime when a fair reference field is 
known. Specifically, we will first assume that the (noise free) data 
is in the image of (9P/9B)g^: 

d = 5PL = {dP /aB) • 5B + e , linear problem (I), 

while for the second problem (the so called Gauss-Newton approx- 
imation) 

d = 5PpL = P — P(Bo) + e , pseudo linear problem (II). 

We investigate the linear problem in this section and the pseudo- 
linear problem in section|4] 



3.1 Linear reconstruction 

Let us illustrate our method on a problem of realistic scales. This 
first simulation is carried on a grid (Nr = 64) with iV^ = 64 
frequencies. The reference field Bo is chosen constant and set to 
1 fiG everywhere for each component, the power spectrum of the 
perturbation field SB has a power law index a = 2 and its RMS 
is 0.01. Data are simulated linearly (see section[3) and are noised 
with a SNR= 20. Figure[T]illustrates the quality of the reconstruc- 
tion. The top panel represents the x and z components along the 
LOS (z direction) or transverse (y direction) for a given pixel. As 



For instance a magnetic loop close to the z axis (where Bx and By ~ 0) 
and its miiTor image by symmetry along the z axis have the same and 
almost zero gradient. 
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Figure 1. Top: input {solid lines) and recovered (dashed lines) x and z components of the field along a LOS (left) and along the y transverse direction (right). 
The y component and the x direction are not plotted since very close to the x component and the y direction. One can see that the z component is better 
reconstructed than the x or y components which is consistant with the variance measurements and the global conditioning of the problem (see section [33) . 
The reconstruction is carried on a A^r = 64 grid with N,, = 64 frequency channel. The data are generated linearly (see section[3) with a SNR= 20. Bottom 
left: maps of |B| for a transverse section after smoothing of the fields. The green images represents the input field while the superposed white contours show 
the recovered one. Bottom right: power spectra of the input field (solid line) and the recovered one (crosses). As expected, the recovered power spectrum is 
damped at higher frequencies because of the regularization. To illustrate this we added the power spectrum of a reconstruction with SNR=200. 



the results for the y component and the x direction are similar to 
the X component and the y direction, they are not plotted. Here, 
the solid lines stand for the input field and the dashed lines for the 
recovered one. It is clear that the two fields are very similar and 
that the z component is the best recovered (see section [33] l. The 
bottom left panel shows a map of for a transverse section af- 
ter smoothing. The smoothing is made by convolving the field with 
a four pixels full-width at half maximum (FWHM) gaussian. The 
green features represent the input field and the reconstructed one 
is shown in the superposed white contours. The bottom right panel 
shows the power spectra of the input field {solid line) and the re- 



covered one {crosses). Finally, figure |2] represents the field lines 
of the input field {top) and the recovered one (bottom). These fig- 
ures show that, if the frequencies are correctly sampled (see section 
13.2b . the linear inverse problem (I) recovers qualitatively well the 
underlying field. The local and global properties of the field can be 
reconstructed provided that the linearization remains valid which 
will be investigated in section|4] 

It is of interest to study the conditioning of the linear prob- 
lem for two reasons (i) to understand the spatial spectral feature 
of the solution; in particular the biases of the eigenvectors of the 
linearized problem which induces anisotropy in the distribution of 
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nential becomes 



2i K B, hm Xi / 




Figure 2. Field lines of the input (left) and the recovered {right) fields for a 
64"^ reconstruction with Ni, = 64 frequencies. The fields correspond to a 
reconstruction with a SNR of 200. 



errors around the solution; (ii) to constrain the best sampling strat- 
egy in order to recover B. Eventually it will also have an impact on 
our ability to carry out the non linear reconstruction. 

The requirements to set up a good conditioning of the global 
inverse problem can be formulated in steps. First a necessary condi- 
tion is to make a proper choice of the (electromagnetic) frequency 
sampling, which can be achieved by looking at a smaller subprob- 
lem on a given LOS; however, this optimal sampling does not war- 
rant a good global conditioning; we therefore investigate the qual- 
ity of the global linear reconstruction by looking at different ele- 
ments of the reconstruction covariance matrix in (spatial) frequency 
space. In particular, we will show that the quality of the reconstruc- 
tion is anisotropic and depends on the components of the field, B, 
which is confirmed by looking at the eigenvectors of the covariance 
matrix for a low dimensional problem. 



3.2 Conditioning of a line of sight and frequency sampling 

One can see easily that in the relation between polarization and 
magnetic field (equation (O), each line of sight is independent of 
the other. The link between them is provided by the solenoidal con- 
dition. In this subsection we will not consider this condition and the 
matrix (9P/9B)bo becomes block-diagonal. Moreover, the three 
components can be separated leading to three different matrices, 
{dP/dB^), (dP/dBy) and (dP/dB^). The field Bo is taken 
constant and its modulus set at 1 ^G. In this case, all blocks are 
the same and the study of the conditioning is reduced to the study 
of three Ni, x A^r matrices with N,^ the number of frequencies and 
A'r the number of pixels in the z direction. 

Numerical investigations show that the conditioning of 
(dP/dBx) depends mainly on the ratio K hric B^/v^ leading to 
the conclusion that the conditioning is dominated by the exponen- 
tial term of equation Q. It follows that {dP/dBy) has the same 
behavior as {dP /dBx) since the exponential terms are the same in 
both equations Q and l|8]l, which is confirmed numerically. 

Recall that since in this section the reference field is chosen 
constant, so is ; therefore the best sampling for the frequencies is 
to have Vn'^ — i^n+i constant, that is a constant step for the squared 
wavelength; hence: = Aq + (n — 1) AA'^ with n — 1, . . . , A'^,^ 
the index of the frequency/wavelength. So that the complex expo- 



(24) 



with m = 1, . . . , A^r the pixel index along the line of sight. The 
value of AA^ must be chosen in such a way that the frequency de- 
pendent complex exponentials are uniformly sampled on the com- 
plex circle. Hence K tic B^ h Nr (n — 1) AA^/c^ must be a multi- 
ple of TV for any n. With L = hNr the maximum probed depth and 
taking the smallest multiple, this yields: 



AA^ 



TT C 



K Uc B, L 



(25) 



With this particular choice, the matrices {dP /dBx) and 
(dP/dBy) take the following form: 



{dP/dB, 



/y) 



IV 



Dtt 



K UcBzL 



X I e "'"e" 



Nr — in 



(26) 



where f3 — 2 K Uc B^ h Xq and C^/y is a different constant in the 

7-1 

X and y directions. If the factor (Aq + (n — l)-!v/{K UcBzL)) ■* 
is set to 1, the matrix is a unitary V andermond matrix and its con- 
ditioning is 1 l lCordovaetaDll99(lh . 

Accounting for this factor impairs the conditioning but it stays 
close to unity. The elements of the last matrix, (dP/dBz) are just 
geometrical series of the elements of {dP /dBx)- Thus, they read: 

T + 3 



1 — exp{—if3{Nr + 1 — m)) exp 



1 w 



KucB^L 
2in 



n-l){Nr + l-m] 



exp 



-i/3exp(-^(n-l))) 



where Cz is yet another constant. At this stage, there is only one 
free parameter left, the first frequency Aq. The conditioning of 
(SP/fBjj/y) being always close to unity, the value of Ao must 
be chosen in order to minimize the conditioning of (dP/dBz). 

Figure^{top patiel) represents the conditioning of (9P/i9B,) 
as a function of Ao for different grid sizes. The curves are very simi- 
lar in shape and the best conditioning is represented by the red dots. 
In the bottom panel the wavelength providing the best conditioning 
for {dP/dBz) is plotted as a function of the grid size. It appears 
that Ao oc \/Nr and the precision on Aq is not really important 
since the minimum of the curves are not really marked. These par- 
ticular choices of Ao give a conditioning of 1.29 for (SP/QB^^/y), 
whatever grid size. 



3.3 Conditioning of Cmap and a posteriori variances 

Let us now investigate the a posteriori variances of different spatial 
frequencies of the reconstructed field. This covariance matrix can 
be written as 



(27) 



where A = (dP/dBo) ■ n with 11 the projector that can- 
cels the divergence (cf. equation J18t ) and C,7^ and = 
/is F~^diag(|fc|")F are the a priori covariance matrices of the 
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Figure 3. Top: conditioning, Cz, of (cJP/cJBz) as a function of Ao for 
different grid sizes. Tlie red dots represent the best conditionings. Bottom: 
Ao giving tlie best conditioning as a function of the grid size, A^r . It appears 
that A() oc VTVr. 

noise and the signal respectivel50. Here we seek Cmap, the Fourier 
transform of Cmap as we want to understand the relative error 
in the amplitude of the spatial modes of B. Because of the po- 
tential high dimensionality of our problem, the covariance matrix, 
Cmap is not computed directly. We chose instead to compute the 
selected values by solving for B the f ollowing equati on with a con- 
jugate gradient method (CGM. IShewchuk 1994 ; Nocedal & Wrighll 

Wmap • B = B,ef. (28) 
Here, Wmap = C^^p and the solution, B, found by the CGM is 

B = CMAP-B,ef. (29) 

The reference field, Brcf, is equal to 1 or ±i for the chosen k fre- 
quency and its opposite — k in order to have a real field, and else- 
where. The elements Bk and B_k of the solution are combinations 
of the covariance of k and — k and the variance of k. It allows us to 
determine the a posteriori variance of the chosen spatial frequency 

^ Throughout this section (unless stated otherwise) we assume that a is 
given by minus the powerspectrum index of the sought magnetic field, and 
choose /Js = 1/P(fc = 1), which corresponds to the minimum variance 
solution. 



k. To check this method, the same variances were also computed 
by the iterative VMLM method. One can check that: 

((Bin - Bout) • (Bin - Bout)^) = Cmap, (30) 

where f denotes conjugate transposition. Bin and Bout stand re- 
spectively for the input field and the reconstructed one in Fourier 
space. As expected, the higher the number of iterations, the closer 
the two estimates of the variance. 

Figure |4]represent the evolution of the a posteriori variance of dif- 
ferent spatial frequencies k for the different components of the field 
in different directions (along a LOS or transverse to it) as a func- 
tion of the SNR. The size of the box is A^r = 16 and the number 
of frequencies is A^^ = 16. Figure |5] shows the evolution of the 
a posteriori variances of the same frequencies as of figure |4l but 
as a function of the spectral index, a, of the sought field (for a 
SNR= 20). As expected, the variance decreases as the index in- 
creases. In Figure|4]the SNR is defined as 

SNR = RMS(data)/crn, (31) 

with CTn Standing for the noise variance. The results for the By and 
Bz fields in the x direction are not plotted because there are exactly 
the same as those in the y direction. First note that the variances, 
af. for the B^ component of the field are much smaller in ampli- 
tude relative to the other components. For the B^; and By fields, at 
low SNR, the Wiener prior is important in the reconstruction, ex- 
plaining the separation of the three curves corresponding to three 
different scales. In Fourier space, Cb ~ fJ-s^ diag(|fcj^") with a 
the spectral index of the power spectrum of the input field. If the 
regularization dominates, Cmap ~ l^l""! which corresponds to 
the values on the figures when the SNR is low. 

For the transverse frequencies {bottom panels), the behaviour 
of the variances is well understood. At low SNR, the Wiener prior 
dominates the reconstruction for the B^: and By components but 
not for the B^ one. Increasing the SNR implies increasing the rel- 
ative weight of the data compared to the prior. So equation i27l 
becomes 

Cmap ~ (A^ ■ C,; ^ ■ A)"\ when SNR ^ oo . (32) 

If we assume a Gaussian white noise, Cn = (Tn I with I the identity 
matrix, equation l l32t becomes 

Cmap ~ cr,^(A^ ■ A.)'', (33) 

so Cmap oc or given equation ( I3U . Cmap oc SNR^^ which 
is the slope of these curves. Finally, note that there is no symmetry 
breaking between the x and y directions and between the x and y 
components of the field or between the sine and cosine modes in 
Cmap- 

Now, consider the x and y components of the field along a LOS 
{top panels). At low SNR, the Wiener prior still dominate, provid- 
ing the same value as in the transverse direction. Then, the variance 
decreases as SNR~^ but reaches a threshold and stagnate. It is clear 
on the figures that there is a symmetry breaking between the x and 
the y components of the field and a separation between the sine and 
cosine modes. At first it may be surprising that the variances reach 
a threshold since the frequencies have been chosen to provide the 
best possible conditioning for (9P/9Bo) along a LOS (see section 
13. 2t . In fact this is a consequence of the solenoidal condition. Re- 
call that for the global inverse problem, the relevant linear model 
is A = (9P/9Bo) ■ n, where II is the projector given by equa- 
tion l llSt . This projector changes the matrix (9P/9Bo) and adds 
off-diagonal terms to the block diagonal matrix considered in the 
previous subsection. In effect, the solenoidal condition degrades 
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Figure 4. A posteriori variance of different spatial frequencies k £ (1, 2, 3) for tlie different components of tlie field in different directions (along a LOS or 
transverse to it) as a function of the SNR. The size of the box is A^r = 16 and the number of frequency is A'^i^ = 16. The top panels con'espond to the variation 
of cr^ for three different values of fc^ while the bottom panels correspond to varying ky . The cosine mode (thick line) and sine mode (dashed line) are both 
shown. All variances decrease with increasing SNR as expected, although at different rate, see the main text. Note the different amplitude in a'^ for the bottom 
right panel which shows that the component of the field is better recovered compared to the other components. This reflects the anisotropy of the model 
A which induces anisotropic reconstruction errors. 



the global conditioning relative to the one LOS problem (but recall 
that without it we have an ill posed problem). In turn this changes 
the eigen structure of Cmap and therefore its projection in Fourier 
space. 

Indeed, let us compute directly the whole matrix Cmap for a 
smaller, more tractable A^r = 8 constant reference magnetic field 
with N,^ = 8 frequencies sampled following the procedure defined 
in section [jjFl Figure |5] shows the global conditioning of the co- 
variance matrix Cmap as a function of the SNR. One can see that 
the mixing of the LOS has a significant effect on conditioning, even 

^ As expected the curves of the variance as a function of the SNR found 
previously are recovered exactly with this direct calculation. 



though the frequencies were chosen optimally. Figure[6]also shows 
that at realistic SNR, the global conditioning remains bounded and 
could be improved, e.g. for the purpose of numerical convergence, 
by artificially increasing the hyperparameter fis- Note finally that 
even though the global conditioning increases with the SNR, the 
variances all decrease, as expected. 



3.4 Eigenspace analysis 

In order to understand the plateau on figure [4] let us also explic- 
itly diagonalize Wmap for the smaller above-described A^r = 8 
problem with a SNR= 20. The corresponding spectrum is plotted 
on figurel?] bottom right panel. The global conditioning of Wmap 
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Figure 5. Same as Figure|4]but as a function of the spectral index, a, of <5B for a SNR= 20. As expected, tlie smootlier tlie expected field, the larger a, the 
smaller the posterior variances. 



is about 10^ (consistently with what was shown on Figure |6] for 
Cmap), but note importantly that there is a cluster of eigenvalues 
followed by a gap. This gap is consistent with the plateau seen on 
figure|4] When increasing the SNR, one expects to filter out less and 
less eigen modes, and therefore to access more and more eigenvec- 
tors (corresponding to decreasing eigenvalues) in the reconstruc- 
tion. However, when reaching the gap, although the SNR increases, 
no more eigenvalues are available for a while. The lower eigenvec- 
tors, encoding informations on higher frequencies, are not within 
reach, and the a posteriori variance of these frequencies stagnate, 
as seen in figure |4] If the SNR increases further, these eigenvalues 
(and therefore their associated eigenvectors) will be sampled, and 
we expect that the al variances will decrease agair0. The modu- 
lus of the first eigenvector (associated to the highest eigenvalue) is 
plotted on the top panels in the x — y (left) and x — z (right) planes. 
It is clear on these figures than the x and y directions are isotropic 
while the z one is anisotropic for this eigenvector. Moreover, the 
component of the power spectra in the bottom left panel show that 
the Bz component clearly differ from the other two components. 

However, all of the main eigenvectors do not behave in the 
same way. Some of them clearly break the symmetry between the x 
and y directions or/and between the x and y components leading to 
the differences in the curves of figure |4] Finally note that the main 
eigenvectors are fairly high frequencies fields. So, the a posteriori 
variances will be smaller for high frequencies than for low ones, 
which is reflected by the top panels of figure |4] 



^ in other words, the plateau seen in the variance per mode in the top panels 

reflects the fact that those modes have non zero contributions from the low 

— 1/2 — 1/2 

signal to noise eigen modes (i.e. eigen modes of • Cb ■ with 

low eigen values, where C^""^ = ■ Cn ""^ ■ A). 




12 3 4 

Log(SNR) 

Figure 6. Global conditioning of the (a posteriori) covariance matrix 
Cmap ^s a function of the SNR. The higher the signal to noise, to more 
difficult the inversion, but the smaller the covariance a posteriori. The 3D 
matrix, A = (SP / 9Bo ) ■ II appears to be more poorly conditioned than its 
ID counterpart even though the sampling in electromagnetic frequency was 
the same as in section 13.21 It remains bounded and within reach of double 
precision calculation. 



4 VALIDITY OF THE LINEAR APPROXIMATION 
4.1 Linear and pseudo linear inversion 

Let us first carry out a linear inversion of the same pertubative field 
5B, with RMS((5B) = 10"^^G, while considering both the linear 
(I) and the pseudo linear (II) data sets (see section[3]l. We work here 
on a A^r = 64 grid, with ~ 64 frequencies, a constant refer- 
ence field of module 1 and SNR=20. Recall that for the linear 
minimum variance solution, the hyperparameter fj,a = 1/P{k = 1) 
(see section l2!2.2t . while for the the pseudo linear data set it may be 
tuned. Figure [8] ?op panel shows the input z component for the in- 
put field (solid line) along a given LOS and the output ones (dotted 
line for the linear data, SPl and dashed line for the pseudo lin- 
ear, 5Ppl) while the bottom panel shows the different power spec- 
tra. As previously, the field recovered from linearized data sets fits 
quite well the input one. The recovered pseudo linear field, though 
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Figure 7. Top panels: maps of the modulus of the field coiTesponding to the first eigenvector of Wmap in th£ G^ft) and x — z (right) plans for a 8^ 
constant reference magnetic field with Ni, = 8 frequencies sampled as explained in section [X2l with a SNR= 20. The first eigenvector appears to be isotropic 
in X and y and anisotropic in the z direction. Bottom left: power spectra of the three components of this eigenvector The anisotropy of the z component is 
clearly visible and in good agreement with the results found in section im cfigurefll and l3.3l (figurel4l. Bottom right: spectrum of the eigenvalues of Wmap. 



somewhat different from the linear one, remains fairly close to the 
original field. The corresponding powers pectra are also shown on 
Figure[8]and confirm that the recovered field in setting (II) is quan- 
titatively redder. 



4.2 Second order residuals 

Let us now study the second order residuals to quantify the do- 
main of validity of the linearization. For this purpose, we subtract 
to the total polarization its zero and first order expansion to obtain 
(P-Po-(9P/9B)3^ oc (5B^) and we divide this quantity by the 
first order term (P — Po oc SH). Figure |9] represents the average 
of this quantity as a function of RMS(5B). Here the perturbation 
consist of a single frequency and single component field. The solid 
lines represent the results obtained with a Bx component along the 
LOS at the lowest mode, while the dashed lines correspond to the 
lowest transverse mode of the Bz component. The dark curi'es rep- 



resent the real part, Q, of the polarization while light ones stand for 
the imaginary part U (see equation (TJ). At very low RMS(5B), 
numerical noise dominate but decreases as the RMS increases. Af- 
ter reaching a minimum, note that the quantity plotted increase as 
RMS((5B) since oc SB'^/SB and thus oc SH. As expected, the 
lower the RMS((5B), the better the linear approximation and the 
better the reconstruction. Note also the significant amplitude dif- 
ference between the B^ and Ji^ components; we interpret this as 
a difference between the second derivatives of the field, which in 
turn, impairs the accuracy of the linearization for the z component. 
This should not be a limitation when carrying the non linear re- 
construction using a method such as VMLM, as the amplitude of 
the subsequent changes in the magnetic field will be scaled by the 
inverse second derivatives. 
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Figure 8. Top: along a LOS for the input field (solid line) and for 
the recovered fields with linear data, SPi^ (dashed line) and pseudo linear 
ones, iPpL (dotted line) (see section [3)- Bottom: power spectra of these 
three fields. Note that the power spectrum of the reconstructed field from 
the pseudo linear data set is steeper. 



4.3 Towards the non linear problem 

Up to now, we have only considered the situation where Bo was 
assumed to be constant. What happens to the conditioning when 
we add spatial frequencies to Bo or/and over Tie? It is easy to see 
that adding transverse frequencies to the x or the y component of 
B will not change the conditioning of a LOS. Indeed, according to 
equations Q and i26\ . only the constants C^/y are modified and 
vary for each LOS, but remain constant along each of them, which 
has no effect on conditioning. On the contrary, if the modulation is 
along a LOS, C^/y is no longer constant, and varies for every pixel 
along a LOS. However, given that the conditioning is dominated by 
the exponential terms in the Vandermond approximation, it doesn't 
change dramatically. Hence the choice of Ao and the sampling fre- 
quency remain the same but the conditioning increases slightly; it 
can reach 3 for {dP/dB^/y) and 40 for (dP/dB,). 
The situation is a priori more dramatic for the z component of 
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Figure 9. Average second order of the polarization divided by the first order 
as a function of RMS((5B). Here RMS(5B) is a single component and 
single mode field. Results are for a the lowest longitudinal mode for the x 
component (solid fines) and the lowest transverse mode for the z component 
(dashed lines). Dark tunes represent the real part Q of the polaiization 
while light ones are for the imaginary part U (see equation Q). 



the field or for the electronic density Tie- Indeed, the addition of a 
transverse modulation has significant consequences, as the value of 
-B, (or/and Tie) in equation ( 1251 ) becomes different for each LOS. 
Therefore, the value of AA^ should in principle be different for 
each LOS to conserve the best conditioning. In practice it is sim- 
plest to take the average of B , (or/and Tie) as a guess. However the 
conditioning per LOS increases signicantly and the quality of the 
reconstruction should be affected. 

However, it appears that the global conditioning of Cmap does not 
change dramatically compared to the constant reference field value, 
whatever the frequency and the amplitude of the added modula- 
tion. The solenoidal condition appears to be very effective. In fact, 
the repetition of the spectral analysis carried in section [34l shows 
that the main difference will be in the gap seen on figure|7] Adding 
modulation on a constant field induces earlier, deeper gaps. At fixed 
SNR, the number of useful eigenvalues for the reconstruction de- 
creases with the modulation. The inversion can still be carried, but 
will be more biased by the lack of resolved eigenmodes. 

As a final illustration, figure [lO] shows an implementation 
of the linear inversion on a more realistic reference field, Bq 
which is e xtracted from a magneto-hydrodynamical simulation 
jKowal & L azarian IQffk. perturbed by a power-law fluctuation 
with a power spectrum of a = 2 and a relative amplitude of 10~^ 
from a virtual data set of SNR=20. Note that for this more realistic 
illustration the electronic density Tie is not constant but extracted 
from the same simulation. Both the shape of the correction and its 
power-spectrum are well recovered for this relative amplitude re- 
flecting that although non constant model and electronic density 
impair the conditionning, reconstructions remain possible. 



5 CONCLUSION AND PERSPECTIVES 

We investigated the problem of reconstructing the three- 
dimensional spatial structure of the magnetic field of a given simu- 
lated patch of our Galaxy, using multi-frequency polarized maps of 
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Figure 10. left top panel: map of a slice (of width 0.047kpc) of input reference magnetic field, Bq; right top panel: map of the same slice but for the known 
electronic density rio; left bottom panel: the input reference magnetic field, the input perturbation and the recovered one along a LOS. The perturbation field 
is a power-law fluctuation with a power spectrum of cs = 2 and a relative amplitude of 10^^ from a virtual data set of SNR=20; right bottom panel: input and 
recovered power spectra of the perturbation field. 



the synchrotron emission at radio wavelengths. 
When starting from a fair approximation of the magnetic field, we 
were able to obtain a good estimate of the underlying field by us- 
ing a linearized version of the inverse problem considered, up to a 
64^ grid size. The spectral analysis of the strictly linear problem 
(with a constant reference field, and the simulated data obtained 
through a linearized model) allowed us to specify the best sam- 
pling strategy in electromagnetic frequency, and predict a spatially 
anisotropic distribution of posterior errors. 

The best sampling strategy is in equal AA^; it follows from the 
shape of (9P/9Bo) along one LOS, which can be approximately 
recast into a unitary Vandermond matrix when this particular sam- 
pling is used. The errors on the reconstructed and By compo- 
nents of the field are shown to be larger than the error on the Bz 
component. This anisotropy can be traced back to the shape of the 



posterior covariance, and ultimately of the linearized model which 
is highly anisotropic, as only the z component of the field induces 
Faraday rotation. 

We considered in turn three more realistic cases: (i) a pseudo 
linear model (linear reconstruction of non-linearly simulated data), 
(ii) a varying reference model Bo, and (iii) a varying reference 
model Bo and a (known) varying electronic density ric- We found 
that for these reconstructions, the global conditioning of the mini- 
mum variance solution remained tractable. Finally, we investigated 
the case where the reference field is given by the outcome of 
a magneto-hydrodynamical simulation, and is perturbed by an 
additional fluctuating component of known power spectrum. We 
showed that even in this case the linear reconstruction quality is 
reasonable. This leads us to claim that a full non-linear reconstruc- 
tion, based on a Gauss-Newton sequence of linear sub-problems of 
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varying reference field, should be achievable. 

Possible extensions of this work, beyond the scope of this 
paper, involve investigating systematically the degeneracies of the 
non-linear inversion. It would be worthwhile to construct specific 
estimators for t he (possibly anisotropic) loc al power spectrum of 
the field (see e.g. lLazarian & Pogosvanl2006l) . Finally, from a mod- 
elling point of view, one of the main limitations of the present 
method is that we had to assume known thermal and relativistic 
electronic densities, in order to obtain a well posed inverse prob- 
lem from synchrotron emission data alone. However, we could in 
principle relax this assumption by adding extra data constrain ing 
the electronic densities (e.g. Hq data, see iHaffner et al. II2OO3') or 
emission measures of pulsars, and attempt a joint reconstruction of 
the magnetic field and the electronic densities. Any prior statisti- 
cal information (e.g. extracted from MHD simulations) of possible 
correlation between B and rio could be used in this context. An- 
other possibility would be to use the extra information given by the 
circular polarization of synchrotron emission (see Appendix |C]l; 
this circular polarization, if negligible in the case of low energy 
sources (like our Galaxy), is measurable in the case of relativis- 
tic radio sources (see e.g. Hone s & Odell 1977), and opens a way 
to constrain the electronic density together with the magnetic field 
structure of the source. 
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APPENDIX A: THE CASE 7 = 3 

For 7 = 3, equation ^ takes a particularly simple expression 



P = A 



J —00 



X exp 

while equation ^ simplifies to 



nAz)(Bi{z) - B;{z) + 2iB^{z)By{z)) 



j\n.B,){z")dz" ]Az, (Al) 



^^il^ = 5^{v-v')2Ahv'^rt,{v'){B^ + iBy){v') 



X exp Y. - ^'){noB.){x' , y' , z")j (A2) 

Note that for this value of 7 the two derivatives with respect to the 
transverse magnetic field are thus related: 



dByiv') dB4r') 



(A3) 



APPENDIX B: SOLENOIDAL FIELDS WITH FIXED 
POWER SPECTRUM. 

The generation of solenoidal (divergence free) fields with fixed 
power spectra up to the Nyquist frequency is a tricky problem. The 
field must obey the three following conditions: 

(i) fixed power spectrum: P(k) oc k^"', 

(ii) free divergence: V-B = 0<4>k-B = 0, 

(iii) reality of the field: Bk ~ B-k* . 

Given conditions (i) and (ii), the field is best generated in Fourier 
space. Since the field is multi periodic and we may write 



B = B±ie±i + B±2ej 



(Bl) 



where ey = k/|k|, exi andex2 form a spherical basis in Fourier 
space, while 13±j, i=l,2 are the projection over that basis of the 
Fourier componant of the field. The vectors exi and e±2 are cho- 
sen in such a way that ekxi/2 ~ — e_kxi/2-The spherical basis is 
direct for k and indirect for — k. In this representation, conditions 
(ii) and (iii) become. 



Bkxi/2 = — S-kxi/2, and BkH 



0. 



(B2) 



So, the first step is to generate two complex fields Bxi and Bx2 
with the sought power spectrum and then apply equation l lB2t . 
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Next, consider the frequencies that have no conjugate, i.e. the fre- 
quency ki = (constant) and ki — Ny (Nyquist frequency) where 
the index i represents the Cartesian coordinates. Let us define Fi 
as the set of these two particular values, i.e. Fi = [Q,Ny], and 
F2 the set of all the other values, i.e. for a vector of dimension A'', 
F2 = [-iN/2-l), -{N/2-2),..., -1,1,. ..N/2-2,N/2~l]. 
When the three components of k belong to Fi, the reality condition 
of the field is merely Tm_B — 0. After putting this imaginary part 
to 0, the field can be projected into the Cartesian basis. 
The difficulty arises when one or two components belong to Fi . For 
example, consider the frequency k = {k^c, ky,kz) with k^ £ Fi, 
ky and k^ £ F2. In this case, condition (iii) become Bk = B_k* 
where k = {kx, —ky, — k^) is the "opposite" of k. The problem 
is that in this case, ekxi/2 7^ ~S-kxi/2 above discussed 

method can no longer apply. Fortunately, the combination of con- 
dition (ii) and Bk = B_^^^ leads to the following set: 



can be reexpressed in terms of the {Q, V, U) "vector" as: 



k B = 0, and Bk, = 0. 



(B3) 



So, the trick is to put the faulty component to and to generate the 
other two as previously but in 2D space. Now, if k2D = {ky,kz), 
we generate B = B±2De±2D, where e||2D = k2D/|k2D| and 
ex2D form a polar basis in Fourier space. As previously, the vectors 
ex2D are chosen in such a way that ekao-Lao = — e_k2D-L2D. In 
this 2D representation, conditions (ii) and (iii) lead to: 



Sk2D-L2D 



^B* 



and B] 



k2Dl|2D 



0. 



(B4) 



Here we have only one degree of freedom left, thus, for these fre- 
quencies, we must generate one complex field -Bx2d with the de- 
sired power spectrum, and then apply equation iB4\ . When ky or 
kz belongs to Fi, a similar procedure applies. 
In the last case, two component belong to Fi. For example, k = 
(kx, ky,kz) with kx G F2, ky and kz G Fi. In this case, condi- 
tion (iii) become Bk = B_k* where k = {—kx,ky,kz) is the 
"opposite" of k. Again, ekxi/2 7^ ^S-kxi/2 ^^'^ combina- 
tion of condition (ii) and Bk ~ B_^^^, leads to equations l lB3b . 
Consequently, the same procedure follows for these frequencies. 
After inverse Fourier transform, one can check that the field is real, 
solenoidal and with the right power spectrum up to the Nyquist fre- 
quency. 



APPENDIX C: CIRCULAR POLARIZATION 

Since the rotating term depends on the density field of thermal elec- 
trons rio in the medium, we cannot separate, with the Faraday ro- 
tation only, rio from B, . One way to tackle this problem is to pick 
up the next coupling term of the Stokes parameters in the (optically 
thin medium, strong rotativity limit) assumption that describes our 
medium. This next term is a factor of conversion between linear and 
circular polarization, that can be considered together with the syn - 
chrotron emissivity of circ ular polar i zation jjones & Odel]||T9T7h . 
Following the notations of ISazono vl j 19691) . we write the transfer 
equation of the polarization tensor 7^^ as follows: 



<ilal3{z) 

dz 



with/„a = 



= Ea0{z) — i{Taa{z)S0-r—SaaTp^{z))Iar{z) , (CI) 

, and Eai3 (z) is an emissivity 



I + Q U + iV 
''-f-\U-iV I-Q 
term. In the assumption of a thin, strongly rotating medium, we can 
retain only the rotating terms (the Hermitian part) of Ta/3. Defining 
y ^ / h q + if 



we can show that the transfer equation 



d_ 

dz 



Q 




Eq 




h 




Q 


V 




Ev 




f 


X 


V 


u 




Eu 




q 




u 



(C2) 



The fact that this differential equation involves multiplication by a 
non-Abelian group element - in SO(3) - prevents us from writing a 
formal solution to the equation in terms of exponentials. However, 
since we are in the end working on a discretized mesh, we can still 
write a formal solution to the discrete problem in terms of (finite) 
sums of (finite) rotations products as we will see below. One impor- 
tant point to notice, linked to the tensor nature of equation ICll is 
the transfoimation law of these "vectors" under rotation of the co- 
ordinate axes in the plane perpendicular to the line of sight. In this 
respect, the vector {h, q, /) behaves the same way as the vector 
(Q, V, U), i.e. the {Q, U)and {h, f) subvectors are rotated by 2-0 
when the coordinate axes are rotated by -0. In the c ase of a homo- 

feneous medium, this allows ISazono 3 l ll969h and Ijones & Odelil 
197% to choose the coordinate axes used to measure Q and U so 
that the V Stokes parameter couples only to U (this is achieved 
when q is set to 0). In this reference frame, the projection of the 
(constant) magnetic field is aligned with the second coordinate axis. 

In the case of a fluctuating magnetic field, such a scheme is not 
possible anymore, and we need to rotate the coupling coefficients 
(best expressed in the reference frame given by the local projec- 
tion of the magnetic field) in a common, constant, reference frame. 
Thus, in an inhomogeneous medium, the equation |C2I in the com- 
mon reference frame takes the form: 

d_ 

dz 

U \ I —-CjQ tillH^^'ip ) J I —lltiLn\^Ztp) J 

(C3) 

where {Q, U, V) are measured in the common reference frame, and 
all other quantities are defined in the frame of the local magnetic 
field. In the applications we will consider in this paper, the rota- 
tion coefficients are dominated by the contribution of cold (thermal) 
electrons^ of the medium. In this context, {h, f) take the following 
form ( iSazonoviil969i) : 



" Q ' 




Bgcos(20) 




h cos(2-0) 




" Q 


V 




Ev 




/ 


X 


V 


u 




-Eq sin(2V') 




—hsm{2tp) 




u 



h = 



qt ric bI 

4 TT^ m'i c-^ 



ql B« 



It is interesting to note that both the frequency dependence, and 
the dependence on the magnetic field are different in the coupling 
terms. We note that IC3I involves the multiplication of the Stokes 
"vector" by an element of a non-Abelian group (SO(3)), which pre- 
cludes finding a formal solution to this differential equation. How- 
ever, the linearity of the equation in the Stokes parameters, allows 
us to write a formal solution in the discretized case in terms of sums 
of products of rotations on the source terms. This equation is very 
similar to the rigid body type equations encountered in mechan- 
ics, with the (major) difference that it is linear. For simplicity, we 
will consider here a first-order discretization of the problem (i.e. we 
consider the different fields to be piecewise constant). The solution 
to the homogeneous Stokes transfer equation can be written as: 



V 

U 



(z) = Yl exp{-AzMi 



V 

U 



(z = 0), 



where A'li is the skew-symmetric matrix corresponding to the vec- 
tor {hi cos{2tl;i), fi, —hi sm{2tpi)). This discretized solution en- 
sures the exact conservation of the polarization degree in the ab- 
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sence of internal sources. It corresponds to the simplest possible 
case of integration of an equation on the SO(3) Lie Group (e.g. 
ICelledoni & OwrenllioOll) . By linearity, we can find the discrete 
solution to the transfer equation with sources l lC3t : 



V 
U 



= I] n exp(-AzM,' 



Egt COs{2i:i) 

Evi 
-EQ,sin{2ipi) 



with ij) = arctan(_B2;/i3j,) 0. This expression generalizes equa- 
tion l|6j and can be used to infer the Frechet derivatives of the po- 
larization field with respect to the magnetic field, as was done in 
section |2] from the integral solution. The source ter ms in the frame 
attached to the local transverse magnetic field read jjones & Odeiil 
Ii923)): 

EQ = Cn,B^v~^ , Ev = ~Dn,B\^Blv~i (C4) 

where n-c is the distribution of high-energy electrons in the medium, 
and C and D are constants that depend on the energy distribution 
of relativistic electrons. Note that in equation l lC4t . n,. is weighted 
differently in the expressions of Eq and Ev, hence we can in prin- 
ciple disentangle B from rir. Here different assumptions can be 
made, namely assuming either that n-c is related to the distribution 
of thermal electrons tIc, or that it is cons tant, or that it is r elated 
to the magnetic field pressure locally (see lBeck et al.l ( |2003|) for a 
discussion of the different assumptions). Another possible path is 
to add external const raints on either no fc oming for instance from 
Hq observations (see lHaffner et al ] |2003i) or from dispersion mea- 



surements of pulsars), or eventually on the relativistic electron dis- 
tribution Tir with diffuse gamma-ray measurements. 



^ Beware that this angle corresponds to 7r — 1/1 in the notations of section|2] 
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